In [1]:
import matplotlib.pyplot as plt
import pandas as pd
import plotly.graph_objs as go
import plotly.express as px
import numpy as np
import os
import sys
import pickle
import seaborn as sns
sns.set_theme(style="whitegrid")
import warnings
warnings.filterwarnings("ignore")

The divergence of the results from reality is due to an overuse of nuclear power (either because EDF has done load following which is not modelled by the programme because nuclear power is too cheap, or because the ENTSOE data gives too much availability) as well as to the non-modelling of the Netherlands and Austria

1. Simulation running or results extraction¶

In [2]:
#Run this for simulation
year=2017 #must be between 2017 and 2019
number_of_sub_techs=5
consumption=pd.read_csv('Input data/'+str(year)+'/'+str(year)+'_MultiNodeAreaConsumption.csv')

if os.path.basename(os.getcwd()) == "Seven_node_Europe":
    os.chdir('../..')
from Main_historical import *

Variables=main_historical(year,number_of_sub_techs)
	 Model creation at 0:00:01.126538
	 Start solving at 0:06:43.273729
	 Solved at 0:14:41.933736
	 Total duration: 0:14:46.614648
	 Model creation at 0:00:01.053781
	 Start solving at 0:07:07.677386
	 Solved at 0:15:06.281685
	 Total duration: 0:15:12.020071
In [3]:
#Run this to load a result file
year=2019 #must be between 2017 and 2019
consumption=pd.read_csv('Input data/'+str(year)+'/'+str(year)+'_MultiNodeAreaConsumption.csv')
with open("Result_"+str(year)+".pickle", 'rb') as f:
    Variables=pickle.load(f)

2. Preliminary calculations¶

In [4]:
states=["Historical"]
exch_list=[Variables["exchange"]]
data=[Variables["energy"]]
data_stor=[Variables["storageOut"]]
data_capa=[Variables["capacity"]]
data_stor_capa=[Variables["Pmax"]]
colors=["black"]
In [5]:
#subtechs list
l=data_capa[0]
l=l.TECHNOLOGIES.unique()
u=[]
n_list=[]
for tech_name in l:
    if tech_name[:4]=="Coal" and tech_name!="Coal":
        n_list.append(tech_name[4:])

3. Results plot¶

a. Capacity¶

In [6]:
# total_consum=round(consum.areaConsumption.sum()/10**6,2)
tech_list=["OldNuke","NewNuke","Coal","Lignite","CCG","TAC","HydroRiver","HydroReservoir",'HydroStorage',"Biomass",
           "WindOnShore","WindOffShore","Solar","Battery"]



df_capa=pd.DataFrame(columns=["TECHNOLOGIES"]+states)
df_capa["TECHNOLOGIES"]=np.array(tech_list)
# print(len(data))
for i in range(len(data)):
    capa=data_capa[i]
    for t in tech_list:
        t_list=[t]
        if len(n_list)>0:
            for j in range(len(n_list)):
                t_list.append(t+n_list[j])

        if t not in ["HydroStorage","Battery"] :
            val=round(capa[capa.TECHNOLOGIES.isin(t_list)].capacity.to_numpy().sum()/10**3,2)
    #         print(val)
            df_capa.loc[df_capa.TECHNOLOGIES==t,states[i]]=val
        else:
            u=data_stor_capa[i]
            val=round(u[u.STOCK_TECHNO==t].Pmax.to_numpy().sum()/10**3,2)
#             print(val)
            df_capa.loc[df_capa.TECHNOLOGIES==t,states[i]]=val
total_capa=df_capa[states[0]].sum()
In [7]:
df_capa
Out[7]:
TECHNOLOGIES Historical
0 OldNuke 98.31
1 NewNuke 0.0
2 Coal 51.5
3 Lignite 21.2
4 CCG 147.11
5 TAC 14.4
6 HydroRiver 29.44
7 HydroReservoir 38.0
8 HydroStorage 38.36
9 Biomass 16.71
10 WindOnShore 114.94
11 WindOffShore 17.32
12 Solar 100.44
13 Battery 0.0
In [8]:
tech_list=["OldNuke","NewNuke","Coal","Lignite","CCG","TAC","HydroRiver","HydroReservoir",'HydroStorage',"Biomass",
           "WindOnShore","WindOffShore","Solar","Battery"]

for area in ["FR","DE","GB","BE","IT","ES","CH"]:
    print("\nCapacities installed for "+area)
    df_capa=pd.DataFrame(columns=["TECHNOLOGIES"]+states)
    df_capa["TECHNOLOGIES"]=np.array(tech_list)
    # print(len(data))
    for i in range(len(data)):
        capa=data_capa[i]
        capa=capa[capa.AREAS==area]
        for t in tech_list:
            t_list=[t]
            if len(n_list)>0:
                for j in range(len(n_list)):
                    t_list.append(t+n_list[j])

            if t not in ["HydroStorage","Battery"] :
                val=round(capa[capa.TECHNOLOGIES.isin(t_list)].capacity.to_numpy().sum()/10**3,2)
        #         print(val)
                df_capa.loc[df_capa.TECHNOLOGIES==t,states[i]]=val
            else:
                u=data_stor_capa[i]
                u=u[u.AREAS==area]
                val=round(u[u.STOCK_TECHNO==t].Pmax.to_numpy().sum()/10**3,2)
    #             print(val)
                df_capa.loc[df_capa.TECHNOLOGIES==t,states[i]]=val
    print(df_capa)
Capacities installed for FR
      TECHNOLOGIES Historical
0          OldNuke      63.13
1          NewNuke        0.0
2             Coal        3.0
3          Lignite        0.0
4              CCG       6.58
5              TAC        0.7
6       HydroRiver      10.96
7   HydroReservoir       8.28
8     HydroStorage       5.02
9          Biomass       1.93
10     WindOnShore      13.61
11    WindOffShore        0.0
12           Solar       8.19
13         Battery        0.0

Capacities installed for DE
      TECHNOLOGIES Historical
0          OldNuke       9.52
1          NewNuke        0.0
2             Coal      25.29
3          Lignite       21.2
4              CCG      17.11
5              TAC      11.42
6       HydroRiver       3.98
7   HydroReservoir        1.3
8     HydroStorage       9.42
9          Biomass       7.75
10     WindOnShore      52.79
11    WindOffShore       6.39
12           Solar       45.3
13         Battery        0.0

Capacities installed for GB
      TECHNOLOGIES Historical
0          OldNuke       9.23
1          NewNuke        0.0
2             Coal       8.79
3          Lignite        0.0
4              CCG       42.2
5              TAC        0.0
6       HydroRiver       1.88
7   HydroReservoir        0.0
8     HydroStorage       2.74
9          Biomass       4.06
10     WindOnShore      13.63
11    WindOffShore       9.38
12           Solar      13.44
13         Battery        0.0

Capacities installed for BE
      TECHNOLOGIES Historical
0          OldNuke       5.94
1          NewNuke        0.0
2             Coal        0.0
3          Lignite        0.0
4              CCG       3.01
5              TAC       2.27
6       HydroRiver       0.17
7   HydroReservoir        0.0
8     HydroStorage       1.31
9          Biomass       0.62
10     WindOnShore       2.25
11    WindOffShore       1.55
12           Solar       3.37
13         Battery        0.0

Capacities installed for IT
      TECHNOLOGIES Historical
0          OldNuke        0.0
1          NewNuke        0.0
2             Coal       4.85
3          Lignite        0.0
4              CCG      47.31
5              TAC        0.0
6       HydroRiver      10.65
7   HydroReservoir       3.86
8     HydroStorage       7.57
9          Biomass       1.54
10     WindOnShore       9.62
11    WindOffShore        0.0
12           Solar      20.87
13         Battery        0.0

Capacities installed for ES
      TECHNOLOGIES Historical
0          OldNuke       7.12
1          NewNuke        0.0
2             Coal       9.56
3          Lignite        0.0
4              CCG      30.38
5              TAC        0.0
6       HydroRiver       1.16
7   HydroReservoir      19.15
8     HydroStorage       5.65
9          Biomass       0.51
10     WindOnShore      22.96
11    WindOffShore        0.0
12           Solar       6.75
13         Battery        0.0

Capacities installed for CH
      TECHNOLOGIES Historical
0          OldNuke       3.37
1          NewNuke        0.0
2             Coal        0.0
3          Lignite        0.0
4              CCG        0.5
5              TAC        0.0
6       HydroRiver       0.63
7   HydroReservoir       5.42
8     HydroStorage       6.64
9          Biomass        0.3
10     WindOnShore       0.08
11    WindOffShore        0.0
12           Solar       2.52
13         Battery        0.0

b. Production¶

In [9]:
country="FR"
p=Variables["capacity"]
p=p[p.AREAS==country]
p
Out[9]:
AREAS TECHNOLOGIES capacity
0 FR Coal6666666666666666 374.625
1 FR curtailment 20000.000
2 FR CCG3333333333333333 823.125
3 FR CCG8333333333333333 823.125
4 FR Lignite3333333333333333 0.000
5 FR Lignite16666666666666666 0.000
6 FR Lignite6666666666666666 0.000
7 FR Biomass3333333333333333 241.375
8 FR TAC8333333333333333 87.875
9 FR TAC 351.500
10 FR Biomass8333333333333333 241.375
11 FR CCG16666666666666666 823.125
12 FR Lignite8333333333333333 0.000
13 FR WindOnShore 13610.000
14 FR Lignite 0.000
15 FR Coal16666666666666666 374.625
16 FR NewNuke 0.000
17 FR Coal3333333333333333 374.625
18 FR WindOffShore 0.000
19 FR Coal8333333333333333 374.625
20 FR CCG6666666666666666 823.125
21 FR TAC16666666666666666 87.875
22 FR OldNuke 63130.000
23 FR HydroReservoir 8279.000
24 FR TAC6666666666666666 87.875
25 FR Biomass16666666666666666 241.375
26 FR Biomass6666666666666666 241.375
27 FR Solar 8188.000
28 FR HydroRiver 10955.000
29 FR TAC3333333333333333 87.875
30 FR Biomass 965.500
31 FR CCG 3292.500
32 FR Coal 1498.500

i. Production at the 7 nodes scale¶

In [10]:
# total_consum=round(consum.areaConsumption.sum()/10**6,2)
tech_list=["OldNuke","NewNuke","Coal","Lignite","CCG","TAC","HydroRiver","HydroReservoir",'HydroStorage',"Biomass",
           "WindOnShore","WindOffShore","Solar","Battery","curtailment"]



df_prod=pd.DataFrame(columns=["TECHNOLOGIES"]+states)
df_prod["TECHNOLOGIES"]=np.array(tech_list)

for i in range(len(data)):
    prod=data[i]
    for t in tech_list:
        t_list=[t]
        if len(n_list)>0:
            for j in range(len(n_list)):
                t_list.append(t+n_list[j])

        if t!="HydroStorage":
            val=round(prod[prod.TECHNOLOGIES.isin(t_list)].energy.to_numpy().sum()/10**6,2)
    #         print(val)
            df_prod.loc[df_prod.TECHNOLOGIES==t,states[i]]=val
        else:
            u=data_stor[i]
            val=round(u[u.STOCK_TECHNO==t].storageOut.to_numpy().sum()/10**6,2)

            df_prod.loc[df_prod.TECHNOLOGIES==t,states[i]]=val
total_prod=df_prod[states[0]].sum()
In [11]:
df_prod_bar=pd.DataFrame(columns=["TECHNOLOGIES","Production (TWh)","Production (%)","Scenario"])
for c in df_prod.columns[1:]:
    u=df_prod[["TECHNOLOGIES",c]]
    u.rename(columns={c:"Production (TWh)"},inplace=True)
    u["Production (%)"]=u["Production (TWh)"]/u["Production (TWh)"].sum()*100
    u["Scenario"]=c
    df_prod_bar=pd.concat([df_prod_bar,u],ignore_index=True)
In [12]:
import plotly.express as px
fig = px.bar(df_prod_bar, x='Scenario', y='Production (%)', color="TECHNOLOGIES")
fig.show()
In [13]:
fig = px.pie(df_prod_bar, values='Production (TWh)', names='TECHNOLOGIES', title='Production per technology (in %)')
fig.show()

ii. Production at country scale¶

In [14]:
# total_consum=round(consum.areaConsumption.sum()/10**6,2)
countries=["FR","DE","BE","GB","IT","CH","ES"]
tech_list=["OldNuke","NewNuke","Biomass","Coal","Lignite","CCG","TAC","HydroRiver","HydroReservoir",'HydroStorage','Battery',
            "WindOnShore","WindOffShore","Solar"]



df_prod_country=pd.DataFrame(columns=["AREAS","TECHNOLOGIES","State","Production (TWh)"])


for i in range(len(states)):
    prod=data[i]
    for c in countries:
        for t in tech_list:
            t_list=[t]
            if len(n_list)>0:
                for j in range(len(n_list)):
                    t_list.append(t+n_list[j])
            if t not in ["HydroStorage","Battery"]:
                val=round(prod[(prod.TECHNOLOGIES.isin(t_list)) & (prod.AREAS==c)].energy.to_numpy().sum()/10**6,2)
                array=np.array([[c,t,states[i],val]])
                df_prod_country=pd.concat([df_prod_country,
                                           pd.DataFrame(array,columns=["AREAS","TECHNOLOGIES","State","Production (TWh)"])])
            else:
                u=data_stor[i]
                val=round(u[(u.STOCK_TECHNO==t) & (u.AREAS==c)].storageOut.to_numpy().sum()/10**6,2)
                array=np.array([[c,t,states[i],val]])
                df_prod_country=pd.concat([df_prod_country,
                                           pd.DataFrame(array,columns=["AREAS","TECHNOLOGIES","State","Production (TWh)"])])

df_prod_country["Production (TWh)"]=df_prod_country["Production (TWh)"].astype('float')
                
# comparison_per_country=df_prod_country.copy()
# for country in countries:
#     total_prod=df_prod[states[0]].sum()
#     total_consum=round(consum[consum.AREAS==country].areaConsumption.sum()/10**6,2)
#     comparison_per_country.loc[comparison_per_country.AREAS==country,"Production (TWh)"]= comparison_per_country.loc[comparison_per_country.AREAS==country,"Production (TWh)"].astype(float)/total_consum
    
# comparison_per_country["Difference with historical"]=np.zeros(len(comparison_per_country))
# for state in states:
#     comparison_per_country.loc[comparison_per_country.State==state,"Difference with "+states[0]]=(comparison_per_country.loc[comparison_per_country.State==state,"Production (TWh)"]-comparison_per_country.loc[comparison_per_country.State=="Historical","Production (TWh)"])*100                
In [15]:
from  matplotlib.ticker import PercentFormatter
g = sns.catplot(x="State", y="Production (TWh)",hue="AREAS",col="TECHNOLOGIES",col_wrap=3,kind="bar",
                data=df_prod_country)
# for ax in g.axes.flat:
#     ax.yaxis.set_major_formatter(PercentFormatter())

plt.show()
In [16]:
country="FR"
p=Variables["energy"]
p=p[p.AREAS==country]
storout=Variables["storageOut"]
#print("HydroStorage",round(storout[storout.AREAS==country].storageOut.sum()/10**6,2),"TWh")
storout=storout[(storout.AREAS==country) & (storout.STOCK_TECHNO!="Battery")]
storout=storout.rename(columns={'STOCK_TECHNO':'TECHNOLOGIES',"storageOut":"energy"})
p=pd.concat([p,storout],ignore_index=True)
tot=0
for t in p.TECHNOLOGIES.unique():
    l=p[p.TECHNOLOGIES==t]
    print(t,round(l.energy.sum()/10**6,2),"TWh")
    tot+=l.energy.sum()
print("Total produced",round(tot/10**6,2),"TWh")
Coal6666666666666666 1.58 TWh
curtailment 0.0 TWh
CCG3333333333333333 2.92 TWh
CCG8333333333333333 0.0 TWh
Lignite3333333333333333 0.0 TWh
Lignite16666666666666666 0.0 TWh
Lignite6666666666666666 0.0 TWh
Biomass3333333333333333 0.0 TWh
TAC8333333333333333 0.0 TWh
TAC 0.0 TWh
Biomass8333333333333333 0.0 TWh
CCG16666666666666666 3.58 TWh
Lignite8333333333333333 0.0 TWh
WindOnShore 32.71 TWh
Lignite 0.0 TWh
Coal16666666666666666 2.16 TWh
NewNuke 0.0 TWh
Coal3333333333333333 2.05 TWh
WindOffShore 0.0 TWh
Coal8333333333333333 1.4 TWh
CCG6666666666666666 0.04 TWh
TAC16666666666666666 0.0 TWh
OldNuke 393.39 TWh
HydroReservoir 17.39 TWh
TAC6666666666666666 0.0 TWh
Biomass16666666666666666 0.0 TWh
Biomass6666666666666666 0.0 TWh
Solar 11.42 TWh
HydroRiver 40.44 TWh
TAC3333333333333333 0.0 TWh
Biomass 0.0 TWh
CCG 1.08 TWh
Coal 4.84 TWh
HydroStorage 2.56 TWh
Total produced 517.56 TWh
In [17]:
tech_list=["OldNuke","NewNuke","Biomass","Coal","CCG","HydroRiver","HydroReservoir",'HydroStorage',
           "WindOnShore","WindOffShore","Solar","curtailment"]
sorted_df=pd.DataFrame(columns=["Date","TECHNOLOGIES","energy"])
for i in range(len(tech_list)):
    d=p[p.TECHNOLOGIES==tech_list[i]].sort_values(by="Date")
    sorted_df=pd.concat([sorted_df,d[["Date","TECHNOLOGIES","energy"]]],ignore_index=True)
fig = px.area(sorted_df, x="Date", y="energy", color="TECHNOLOGIES",line_group="TECHNOLOGIES",title=country+" simulations")
consum=consumption[consumption.AREAS==country]
fig.add_scatter(x=consum.Date, y=consum.areaConsumption, name="Consumption",mode="lines",line_color="black")
fig.show()

c. HydroStorage operations¶

In [18]:
stocklevels=[Variables["stockLevel"]]
for area in stocklevels[0].AREAS.unique():
    fig, ax = plt.subplots(1,figsize=(20,10))
    fig.suptitle(area+" Hydrostorage stock level")
    v=stocklevels[0]
    u=v[v.AREAS==area]
    u=u[u.STOCK_TECHNO=="HydroStorage"]
    ax.plot(u.Date,u.stockLevel,"o",color=colors[0],label=states[0])
    ax.legend()
    ax.set_ylabel("Storage (MWh)")
    plt.xlabel("Date")

d. Exchanges¶

In [19]:
c1,c2="FR","ES"

plt.figure(figsize=(10,5))

fig=plt.figure(figsize=(20,5))
plt.title(c1+'-'+c2+" exchanges")
    
u=exch_list[0]
v=u[(u.AREAS==c1)&(u.AREAS1==c2)]
v.TIMESTAMP=pd.to_datetime(v.TIMESTAMP)
plt.scatter(v.TIMESTAMP,v.exchange,label=states[0],s=20,alpha=0.3,color=colors[0])
plt.legend()
plt.ylabel('Exchanges (MW)')
# plt.subplots_adjust(left=None, bottom=0, right=None, top=1, wspace=0.2, hspace=None)
plt.xlabel("Date")
Out[19]:
Text(0.5, 0, 'Date')
<Figure size 720x360 with 0 Axes>
In [20]:
ext_fr=["DE","BE","GB","IT","CH","ES"]
exchanges_fr=pd.DataFrame(columns=["state",'countries',"Exchanges (TWh)"])
other_exchanges=pd.DataFrame(columns=["state",'countries',"Exchanges (TWh)"])

for i in range(len(exch_list)):
    s=0
    exchange=exch_list[i]
    for c in ["DE","BE","GB","IT","CH","ES"]:
        value=round(exchange[(exchange.AREAS=="FR") & (exchange.AREAS1==c)].exchange.sum()/10**6,2)
        s+=value
        u=pd.DataFrame(np.array([[states[i],c,
                                 value]]),
                       columns=["state",'countries',"Exchanges (TWh)"])
        exchanges_fr=pd.concat([exchanges_fr,u],ignore_index=True)

    u=pd.DataFrame(np.array([[states[i],"Total",
                                 s]]),
                       columns=["state",'countries',"Exchanges (TWh)"])
    exchanges_fr=pd.concat([exchanges_fr,u],ignore_index=True)

    for a,b in [["CH","DE"],["CH","IT"],["DE","IT"]]:
        value=round(exchange[(exchange.AREAS==a) & (exchange.AREAS1==b)].exchange.sum()/10**6,2)
        u=pd.DataFrame(np.array([[states[i],a+"-"+b,
                                 value]]),
                       columns=["state",'countries',"Exchanges (TWh)"])
        other_exchanges=pd.concat([other_exchanges,u],ignore_index=True)
        
exchanges_fr["Exchanges (TWh)"]=exchanges_fr["Exchanges (TWh)"].astype(float)
other_exchanges["Exchanges (TWh)"]=other_exchanges["Exchanges (TWh)"].astype(float)
In [21]:
ax = sns.barplot(x="countries", y="Exchanges (TWh)", hue="state",data=exchanges_fr)
ax.title.set_text("France exchanges")
In [22]:
ax = sns.barplot(x="countries", y="Exchanges (TWh)", hue="state",data=other_exchanges)
ax.title.set_text("Other exchanges")

e. Cost¶

In [23]:
capacityCost=np.array([Variables['capacityCosts'].capacityCosts.sum()/10**9])
In [24]:
energyCost=np.array([Variables['energyCosts'].energyCosts.sum()/10**9])
In [25]:
ax=sns.barplot(x=["Scenario"],y=capacityCost)
plt.ylabel("B€/year")
plt.title("Annualized capacity cost")
ax.bar_label(ax.containers[0])
Out[25]:
[Text(0, 0, '102.442')]
In [26]:
ax=sns.barplot(x=["Scenario"],y=energyCost)
plt.ylabel("B€/year")
plt.title("Annualized energy cost")
ax.bar_label(ax.containers[0])
Out[26]:
[Text(0, 0, '31.3012')]
In [27]:
ax=sns.barplot(x=["Scenario"],y=energyCost+capacityCost)
plt.ylabel("B€/year")
plt.title("Annualized total cost")
ax.bar_label(ax.containers[0])
Out[27]:
[Text(0, 0, '133.743')]
In [28]:
tot_cost=energyCost+capacityCost
ax=sns.barplot(x=["Scenario"],y=tot_cost*10**9/np.array([total_prod])*10**-6)
plt.ylabel("€/MWh")
plt.title("Cost of electricity")
ax.bar_label(ax.containers[0])
Out[28]:
[Text(0, 0, '68.0256')]